

2,e- /7? y 


NASA-CR-1 78038 
19860009595 

ICASE 


NASA Contractor Report 178038 
ICASE REPORT NO. 85-59 


ACCURACY OF SCHEMES FOR THE EULER EQUATIONS 
WITH NON-UNIFORM MESHES 


E. Turkel 
S. Yaniv 
U . Landau 




i 


Contract Nos. NAS1-17070 and NAS1-18107 
December 1985 


INSTITUTE FOR COMPUTER APPLICATIONS IN SCIENCE AND ENGINEERING 
NASA Langley Research Center, Hampton, Virginia 23665 

Operated by the Universities Space Research Association 


NASA 

National Aeronautics and 
Space Administration 

Langley Research Center 

Hampton, Virginia 23665 




ACCURACY OF SCHEMES FOR THE EULER EQUATIONS WITH NON-UNIFORM MESHES 


E. Turkel 

Institute for Computer Applications in Science and Engineering 

and 

Tel-Aviv University 


S. Yaniv and U. Landau 
IAI, Israel 


ABSTRACT 

The effect of non-uniform grids on the solution of the Euler equations is 
analyzed. We consider a Runge-Kutta type scheme based on a finite volume 
formulation. We show that for arbitrary grids the scheme can be inconsistent 
even though it is second-order accurate for uniform grids. An improvement is 
suggested which leads to at least first-order accuracy for general grids. 
Test cases are presented in both two- and three-space dimensions. 
Applications to finite difference and implicit algorithms are also given. 
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1. INTRODUCTION 


In recent years much progress has been made in the solution of the steady- 
state Euler equations in both two and three dimensions. For complex shapes 
these calculations are usually based on a body-fitted curvilinear grid. For 
general three-dimensional bodies it is very difficult to construct a body- 
fitted coordinate system, see e.g., [3]. In particular, the grid system used 
frequently is not smooth. In some cases there may exist discontinuities in 
the gradients of the grid while in other cases the gradients vary sharply but 
not discontinuously. 



Figure 1 
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In Figure 1 we present a typical grid constructed by FL057 for a wing-alone 
configuration. This is an expanded view of the trailing edge for one slice in 
the spanwise direction. We see that even in this simple case there are large 
variations between neighboring cells. 


2. SPACE DISCRETIZATION 

We consider the Euler equations for an inviscid fluid written in 
conservation form. For simplicity of representation we shall consider two- 
dimensional flow, however, all the results generalize to three dimensions in a 
straightforward manner. We thus consider 

w t + f x + g y = 0. (1) 

where w = (p, pu, pv, E) 1 ". Rewriting (1) in integral form, we get 

■jj— ff w dv + J F*nds = 0 (2) 

dt V S 

and n is the outward normal. Letting w be the cell- 
w, we can rewrite (2) as 

v -^+J + J + / + J (fdy - gdx) = 0 (3) 

AB BC CD DA 

for the zone shown in Figure 2. At this stage, (3) is still exact. In order 
to solve (3), we now replace the integrals in (3) by an integration rule. 


where F = (fjg)* 1 
averaged values of 
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Figure 2 


Midpoint Integration Rule 

Using the formulation presented in [2] we assume that the metrics are 
defined at the mesh nodes while the dependent variables are cell averages or 
else are located in the center of the cell. Replacing the integrals in (3) by 
the midpoint rule we get 


9w 

v + q e + q f + % + q h - 0 <4) 


where Q = fAy - gAx. We examine this more carefully by considering the 
point G = (i+ V2 » j)> then 
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Q G = f G (y C " y B } “ g G (x C ~ X B } ’ (5 > 

We are then left with the job of evaluating fg, gg. In the original 
formulation [2], w was calculated at G by 



and then fg = f(wg), gg = g(wg). Thus, both cells contribute equally to the 
value of w along the face independent of the geometry of the grid. Thus, we 
do not account for non-uniform volumes and stretchings and deformations in the 
various directions. Another possibility is to average f and g directly. 

f G = 2 ^ f i+l,j + f i,j j * (6a) 

These two procedures are equivalent to within second-order accuracy. There 
are some minor differences with respect to shock resolution. In order to 
check the accuracy of (6) we first consider a one-dimensional example. 

The one-dimensional equivalent of (4) - (5) (see Figure 3) is 



Figure 3 
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Using the averaging on f given by (6a) we get 



dt 


j+1 


- f 


2h . 
J 


J-i = 0. 


( 8 ) 


A straightforward Taylor-series expansion shows that (8) is inconsistent 
unless 


h 1+l " 2h 1 + h i-l _ 


h. 

J 


= 0(h) 


(9) 


where h = max(hj). Let r^ = h j+]/ h j* In [6] we define a grid as algebraic 
if r j = 1 + 0(hP) and exponential if rj + j/rj = 1 + 0(h). Thus, (8) is 
consistent only if the grid is algebraic, and is second-order accurate only if 
we also have p = 2. 

We next replace (6) by 


W 


h i Vi + h i+i “1 

h j + h j+l 


(10) 


and (8) by 


dw f(w.. i, ) - f(w. 1/ ) 

J. + — Hill— Hlh—.a. ai) 

dt ^ 

This scheme is now first-order accurate and furthermore, it is second-order 
accurate for all algebraic grids. In [6] it is also shown how to construct a 
second-order finite difference version for exponentially stretched meshes. 

Until now we have discussed the accuracy of the flux differences. In the 
Runge-Kutta schemes [2,5], it is also necessary to add an artificial 
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viscosity. This consists of two parts: one is an approximation to a second 
d^ff er 6nce with a nonlinear coefficient that is of higher order in smooth 
regions. This essentially enforces an entropy-like condition and ensures the 
correct shock conditions. The second portion of the viscosity is an approxi- 
mation to a fourth difference and is a linear term. This term is needed in 
smooth regions to prevent decoupling the odd and even points and to suppress 
the highest frequencies [5]. These viscosities are now also used in implicit 
A.D.I. algorithms [4] and also in Lax— Wendroff type methods [1]. Both these 
viscosities add third-order terms to the truncation error when the grid is 
sufficiently smooth. However, for rapidly varying grids these viscosities can 
reduce the order of the method. In particular, we have found the lift 
coefficient to be sensitive to the coefficients of these viscosity terms. 


Implementation — Two-Dimensions 

The formula given by (10) is a one-dimensional formula. In order to 
generalize to two dimensions we use this formula in each direction. Using the 
original mapping we calculate the position of the center of the cell and also 
the center of each face. Let (£»ri) be a plane for which the grid is uniform 
and (x,y) be the physical plane. We assume that the mesh is given by a 
mapping T, (x^y^) = TCg^n^). Then the center of a cell is defined by the 
image of (S ±+ l ^ , n.- + y ) while the center of the faces are given by the 
image of (? i± i^ , n^) and rij ± i^). We note that using the explicit 
mapping of [3] the mapping must be adjusted so that the center of the outer- 
most cell is defined and is not at infinity. Having found these positions we 
calculate the distance from the center of the zone to the center of each face. 
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This gives Che distance £ and ra shown in Figure 4. Formula (10) is then 
replaced by 


f G = 


£f . . , . + mf . . 
1+1 >3 1.3 

£ + m 


( 12 ) 


= a f... . + (l-a)f. 
i+l,J i.y 


with a = £/(£+m). 




Figure 4 


We see from (12) that only the ratio £/(£+m) need be known for each zone. 
Hence, once the cell centers and surface centers are found, then only one 
value need be stored per cell per direction. Thus, in three dimensions we 
need to store three three-dimensional vectors and then the positions of the 
cell and face centers can be discarded. 

In the codes used, FL052, FL053 , FL057 , FL059 the metrics are given 

explicitly by analytic formulae. Hence, it is easy to find the cell and face 
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centers. For other mesh generation techniques it is not as easy to explicitly 
find these other positions. One possibility is to construct a grid with twice 
as many grid lines. Then the nodal positions are found by keeping every other 
point. The other grid points on the fine grid give the cell and face centers 
on the final grid. Since we do not need to store all this information, this 
is not a large cost. A more approximate procedure is to calculate the cell 
and face centers by simple averaging from the nodal positions. When the grid 
is highly stretched in the outer field one must be careful calculating the 
averages in the far field. 

We note that the formulas given by (12) are one-dimensional formulas given 
along each coordinate line given by the mapping and not along x and y 
coordinate lines. Hence, this improvement accounts for non-uniform 
stretchings but does not necessarily account for shearing effects between 
cells. Furthermore, if we wish a constant solution to be preserved in the far 
field, it is necessary to average w by (12) and then calculate f. 
Averaging f will destroy constant solutions since the metrics will be 
averaged differently in each direction. 


Trapezoidal Integration Rule 

In the previous sections we replaced each line integral by a midpoint 
rule. Thus, (see Figure 2) 

C 

/ fdy- gdx * (fAy - gAx),,. (13) 

B G 


For a general mesh this formula is only first- order accurate since G is 
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defined as the image of (£ . i, , n.) and is not necessarily the center of 

i+ 72 3 

the line BC. Choosing G to be the center of the line BC would cause trouble 
with the interpolation formulas, e.g., (12). An alternative is to replace the 
midpoint rule (13) by the trapezoidal rule 


C 1 
/ fdy - gdx - j 


Ay(f(B) + f(C)) - Ax(g(B) + g(C)) 


(14) 


This formula is now second-order accurate for non-uniform meshes. In order to 
evaluate f at node, e.g., B, we use bilinear interpolation. This interpola- 
tion Is done in a standard square. Thus, we know the value of w at the cell 

centers (x i} j), (x i+1>j , y i+ i >y ), yi,j+l ) and (x i+l,j+l» 

y^ + ^ j+i)» We now use an isoparametric bilinear mapping to map the 
quadrilateral given by these points into the unit square in (£,ri) space. 
Hence, 


x = Aj + A 2 £ + Ag n + £n (15a) 

y = B^ + B^ K + Bg n + B^ 5 n. (15b) 

We then assume that w satisfies the same mapping, so 

w = c i + c 2 5 + c 3 n + c 4 5n* (15c) 

Knowing x, y, w at the four cell centers corresponding to g,q = 0,1 we can 
evaluate the coefficients A^, B i} We next calculate the value of £ and 

n that gives the values of x and y for node B. Given that value of 
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5 and Ti we calculate w at the node B by (15c). As is standard in finite 
element theory these mappings preserve the continuity of x, y and w between 
zones and the interpolation is first-order accurate. 

This method is more accurate than the midpoint rule previously discussed 
but is also more time consuming. Further, the bilinear transformation needs 
to be modified near boundaries. 

In the case of a uniform grid the midpoint rule and the trapezoidal rule 
give rise to different schemes for advancing w in time. For simplicity of 
notation, we assume that the fluxes are averaged rather than w. For a 
uniform grid the midpoint rule reduces to 


9w ij + f i+l,j f i-l,j g i,j+l gj,,i-l _ 

3t 2Ax 2 Ay 


(16) 


For the uniform grid the trapezoidal rule reduces to 


9w 


1 


3t 8Ax 


f i+l,j+l + 2f i+l,j + f i+l,j-l 


- (f . , ...+21.. . + f . . . .) 
v i-l,j+l i-l,j i-1 , j-1 J 


+ 8Ay [ s i+l,j+l + 2g i,j+l + s i-l , j+1 


" (g i+l,j-l + 2g i,j-l + g i-l,j-J 


= 0 . 


(17) 


We now analyze the stability of both (16) and (17). Linearizing the equation 
and freezing coefficients, we obtain 
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V7 + Aw + Bw =0. 
t x y 

Using a Leapfrog or Runge-Kutta scheme the stability depends on 
values of the matrix 


D = ccA + BB 


where 


a = sin 6, g = sin <p 

. a /l + cos a\ . I l + cos s\ 

a = sin 0l ^ T ) > 3 = sin <f>( 2 / 


We consider a general coordinate system (£,n) and denote the 
coordinates by (x,y). To simplify the analysis we symmetrize D 
find (see [5]) that 



where c is the speed of sound and 


a = ay - By_, b = Bx c 

n 5 5 


ax , w = aq + Br 
n 


q = y u - x v and r = x_ v - y_ u. 

n n 5 5 


(18) 
the eigen- 

(19) 

for (16), 

(for (17). 

Cartesian 
. We then 

( 20 ) 


with 


( 21 ) 
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Hence, the eigenvalues of D are given by 


d = W, W ± 


2 2 
a + b c 


( 22 ) 


It follows that the time step restriction for (16) is of the form 


At K 

V | q | + | r | + L x c 


(23) 


where K is a scheme-dependent constant, and 


Li = + x 


2 2 2 L 2 7 2 2 

n + y C * y n + V x ( + y (~ x n’ 


y n )2 + £U C V y s y n ) ' 


For an orthogonal mesh, L reduces to 


t o ( 2 . 2 2 2^ 

L ! = 2.max(x § + y ? , x^ + y^J. 


The time restriction for the trapezoidal rule (17) is of the form 


and 


At 



K 



\T < 

111 

+ 

r max( q 

LJ 



(24) 


+ L 2 c 


27 

L 2 " 32 V 


Thus, the time restriction for the trapezoidal rule is slightly less stringent 
than that for the midpoint rule. The coupling of modes between mesh points 
also differs between the midpoint and trapezoidal rules. 
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3. COMPUTATIONAL RESULTS 

We first consider the two-dimensional code FL052. This uses a Runge-Kutta 
algorithm for an 0 mesh with acceleration by local time step, enthalpy 
damping, and residual smoothing [2,5] but without any multigrid techniques. 
All calculations were done on a coarse 64x16 mesh. As the mesh is refined the 
accuracy of all the methods gets better. 

The first case that we consider is flow about a NACA 0012 airfoil with 
M, = 0.3 and a = 10 . The resultant flow is subsonic everywhere and so the 
solution could be computed using the full potential equation. The lift 
coefficient is then found to be about 1.27. In Table 1 we present some 
results for this case. In the first column we indicate which case is being 
presented. In the second column we present the coefficient of the fourth- 
order viscosity. In the third and fourth columns we present and Cq 

respectively. In the last column we give the RMS of the density residual 
after 1000 time steps. All calculations were done with a four-step Runge- 
Kutta method with a x = 1/4, a 2 = 1/3, Oj = 1/2, a 4 = 1. This gives fourth- 
order accuracy in time when the coefficients of the differential equation are 
independent of time. Since the Euler equations are nonlinear there is only 
second-order time accuracy. In any case, the time accuracy is not of 
importance since we are considering only steady flows. 

Since the flow is subsonic the second-order viscosity does not play an 
important role. We see from Table 1 that accounting for the non-uniform mesh 
increases the accuracy of both the lift and drag coefficients. We also see 
that for this coarse mesh that the value of VIS4 also greatly affects the 
accuracy. All these calculations were done with a distribution of points in 
the normal direction as given in the original version of FL052. Other runs 
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not shown were done with an exponential-like stretching in the normal 
direction. This provides a smoother mesh and less severe stretching in the 
far field. As expected the use of (10) does not improve things as much when 
the mesh is smoother. In fact, for a uniform mesh (6) and (10) coincide. In 
all cases computed the use of (10) improved the accuracy of the solution. 

The next case that we consider is flow about a NACA 0012 airfoil with M = 
0.85 and a = 1°. In this case there is a strong shock. Comparisons with 
other codes indicates that the C L ~ 0.36 - 0.40 and C D ~ 0.55. As before, 
the lift coefficient is underpredicted by the coarse mesh. Increasing VIS4 
increases C L and also slightly increases Cp. As before we find that the 
use of the weighting formula (10) increases the accuracy of both and 
C D . In this case, we see that the weighted formula also substantially 
increases the convergence rate in some cases. As before, using a smoother 
stretching in the normal direction increases the accuracy of all the cases, 
though using (10) still marginally increases the accuracy. In this case, we 
again see that the value of VIS4 affects accuracy while other computations 
show that VIS2 can also have a strong effect on the accuracy of the code. 

As the final two-dimensional inviscid calculation we consider flow about 
an RAE 2822 airfoil with M = 0.75 and a = 3°. This is also a transonic 

CO 

case which needs good shock resolution. Runs with finer meshes indicate that 
C L - 1.10 and C D ~ .042. In this case the C L is again underpredicted by 
10-20% and the use of the weighting formula (10) again increases the accuracy 
as does increasing VIS4. It is to be noted that the fine mesh calculations 
are still done with the Euler equations and we are not considering viscous 
effects. In this case the use of a smoother exponential-like stretching in 
the normal direction again dramatically increases the accuracy but at the 
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expense of a slower convergence rate. 

We next tested the weighting formulas in three dimensions using FL057 for 
a wing-alone configuration. We found that the code still converged for 
extreme cases as M = 0.60 and a = 24®. The code was run with a coarse 

mesh of 64x10x10 and a finer (though still coarse) mesh of 96x16x20. The 
grid generation used the standard geometry routines contained in FL057. 
Typical results are shown in Table 4. The use of small value for VIS4, less 
than 0.5, introduced oscillations in the solution and sometimes the method 
would not even converge. The use of (10) in each- coordinate direction allowed 
for faster convergence and fewer oscillations in the converged solution. As 
VIS4 was increased to 1.0 the oscillations became less important. Also, as 
the mesh was refined the oscillations diminished. For the case M = 0.60, 
a = 24 and the 96x16x20 mesh, the standard method would not converge 

unless we set VIS2 = 1.0. Using the weighting formula (10) the code coverged 
with VIS2 = 0.5. Recalculating the midpoints of the cells and cell surfaces 
at each time step adds about 30/£ to the running time. If the distance ratios 
in each direction are stored then the increase in running time is negligible 
but three additional three-dimensional vectors need to be stored. 


4. CONCLUSIONS 

Using a Taylor-series expansion one can show that the standard weighting 
of the fluxes (6) at interfaces for a finite volume scheme can be even 
inconsistent for general meshes. A new weighting formula (10) guarantees that 
the scheme is at least first-order accurate for all meshes. As expected, 
computations demonstrate that for smooth meshes the differences between (6) 
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and (10) are minimal. When the stretchings are more severe then (10) can 
offer significant advantages, especially on coarser meshes. In addition, in 
three dimensions we found that the weighting (10) can also improve the 
convergence properties of the scheme. We also found that the constants 
multiplying the artificial viscosity can have a large effect on the accuracy 
of the solution, especially on the lift coefficient. Weighting formulas for 
nodal schemes and MacCormack schemes are given in [6]. 
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TABLE Is 

Flow about 
0 mesh with 

NACA 0012 with 
64 x 16 grid. 

M = 0.3 , a = 
00 

10°. 

Weighting 

VIS4 

C L 

C D 

Residual 

(6) 

0.2 

1.227 

.008 

9 x 10" 9 


0.5 

1.240 

.007 

6 x 10~ 9 

(10) 

0.2 

1.264 

.002 

2 x 10" 8 


0.5 

1.288 

-.001 

1 x 10" 8 

finer mesh 


-1.274 

0 



TABLE 2: 

Flow about NACA 
0 mesh with 64 x 

0012 with 
16 grid. 

H = 0.85, a - 1° 
00 7 


Weighting 

VIS4 

C L 

C D 

Residual 

(6) 

0.2 

.2631 

.0487 

1 X 10“ 9 


0.5 

.2738 

.0499 

1 x 10" 5 

(10) 

0.2 

.2772 

.0492 

1 x 10“ 10 


0.5 

.3006 

.0513 

2 x 10" 8 

finer mesh 

rsA 

.36-. 40 

-.054 
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TABLE 3: Flow about RAE 2822 with M = 0.75, a = 3°. 

OO 


Weighting VIS4 C L C D Residual 


(6) 

0.2 

.9586 

.0395 

H 

l 

o 

r-4 

X 

-H 


0.5 

.9910 

.0407 

5 x 10" 12 

(10) 

0.2 

1.003 

.0414 

3 x 10 -11 


0.5 

1.059 

.0431 

5 x 10" 12 

; iner mesh 


-1.10 

-.042 



TABLE 

4: 

Th r e e-D ime ns i o na 1 

Case with M 

OO 

= 0.60, a = 

24°. 

Weighting 

VIS4 

C L 

C D 

No. 

Iterations 

Residual 

(6) 

0.25 

1.0248 

.4170 

300 

8 x 10 -5 


1.00 

1.1359 

.4611 

150 

7 x 10 -6 


(10) 0.25 1.0155 


.4132 300 


4 x 10 
7 x 10" 6 


1.00 


1.1265 


4571 


150 






2. Government Accession No. 


3. Recipient'^ Catalog No. 


i *'**' No NASA CR- 178038 
ICASE Report No. 85-59 


SIAM J. Sci. Stat. J. 


16. Abstract 

The effect of non-uniform grids on the solution of the Euler equations is 
analyzed. We consider a Runge— Kutta type scheme based on a finite volume 
formulation. We show that for arbitrary grids the scheme can be inconsistent 
even though it is second-order accurate for uniform grids. An improvement is 
suggested which leads to at least first-order accuracy for general grids. 
Test cases are presented in both two- and three-space dimensions. 
Applications to finite difference and implicit algorithms are also given. 


For sale by the National Technical Information Service, Springfield, Virginia 22161 


18. Distribution Statement 

64 - Numerical Analysis 
02 - Aerodynamics 

Unclassified - Unlimited 


21. No. of Pages 

22. Price 

21 

A0 2 


17. Key Words (Suggested by Author(s)) 

accuracy 
Euler equations 
Runge-Kutta scheme 
non-uniform grids 

20. Security Classif. (of this 

Unclassified 


19. Security Oassif. (of this report) 

Unclassified 


5. Report Date 

December 1985 

6. Performing Organization Code 

8. Performing Organization Report No. 




10. Work Unit No. 

11. Contract or Grant No. 

■NAS 1=17070; -NAS 1-18107 

13. Type of Report and Period Covered 

Contractor Report 

14. Sponsoring Agency Code 


4 Title and Subtitle 

ACCURACY OF SCHEMES FOR THE EULER EQUATIONS 
WITH NON-UNIFORM MESHES 

7. Author(i) 

E. Turkel, S. Yaniv and U. Landau 

9. Performing Organization Name and Addreu 

Institute for Computer Applications in Science 
and Engineering 

Mail Stop 132C, NASA Langley Research Center 


ifjinnw<Ti 


12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 
Washington, D.C. 20546 

15. Supplementary Notes 

Langley Technical Monitor: Submitted 

C. South Jr. Comput. 

Final Report 



















